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Abstract 

We consider the problem of learning sparsely used dictionaries with an arbitrary square 
dictionary and a random, sparse coefficient matrix. We prove that 0{n log n) samples are 
sufficient to uniquely determine the coefficient matrix. Based on this proof, we design 
a polynomial-time algorithm, called Exact Recovery of Sparsely-Used Dictionaries (ER- 
SpUD) , and prove that it probably recovers the dictionary and coefficient matrix when the 
coefficient matrix is sufficiently sparse. Simulation results show that ER-SpUD reveals the 
true dictionary as well as the coefficients with probability higher than many state-of-the-art 
algorithms. 

Keywords: Dictionary learning, matrix decomposition, matrix sparsification. 
1. Introduction 

In the Sparsely-Used Dictionary Learning Problem, one is given a matrix Y E M. nxp and 
asked to find a pair of matrices A 6 ]& nxm an d x G M. mxp so that \\Y — AX\\ is small and 
so that X is sparse - X has only a few nonzero elements. We examine solutions to this 
problem in which A is a basis, so m = n, and without the presence of noise, in which case 
we insist Y = AX. Variants of this problem arise in different contexts in machine learning, 
signal processing, and even computational neuroscience. We list two prominent examples: 

• Dictionary learning [17; 13]: Here, the goal is to find a basis A that most compactly 
represents a given set of sample data. Techniques based on learned dictionaries have 
performed quite well in a number of applications in signal and image processing [3; 
19; 22]. 

• Blind source separation [24]: Here, the rows of X are considered the emissions of var- 
ious sources over time. The sources are linearly mixed by A (instantaneous mixing). 
Sparse component analysis [24; 9] is the problem of using the prior information that 
the sources are sparse in some domain to unmix Y and obtain (A, X). 

These applications raise several basic questions. First, when is the problem well-posed? 
More precisely, suppose that Y is indeed the product of some unknown dictionary A and 
sparse coefficient matrix X. Is it possible to identify A and X, up to scaling and permu- 
tation. If we assume that the rows of X are sampled from independent random sources, 
classical, general results in the literature on Independent Component Analysis imply that 
the problem is solvable in the large sample limit [4] . If we instead assume that the columns 
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of X each have at most k nonzero entries, and that for each possible pattern of nonzeros, we 
have observed k + 1 nondegenerate samples yj, the problem is again well-posed [14; 9]. This 
suggests a sample requirement of p > (k + We ask: is this large number necessary? 

Or could it be that the desired factorization is unique 1 even with more realistic sample 
sizes? 

Second, suppose that we know that the problem is well-posed. Can it be solved effi- 
ciently? This question has been vigorously investigated by many authors, starting from 
seminal work of Olshausen and Field [17], and continuing with the development of alternat- 
ing directions methods such as the Method of Optimal Directions (MOD) [5], K-SVD [1], 
and more recent, scalable variants [15]. This dominant approach to dictionary learning ex- 
ploits the fact that the constraint Y = AX is bilinear. Because the problem is nonconvex, 
spurious local minima are a concern in practice, and even in the cases where the algorithms 
perform well empirically, providing global theoretical guarantees would be a daunting task. 
Even the local properties of the problem have only recently begun to be studied carefully. 
For example, [11; 8] have shown that under certain natural random models for X, the 
desired solution will be a local minimum of the objective function with high probability. 
However, these results do not guarantee correct recovery by any efficient algorithm. 

In this work, we contribute to the understanding of both of these questions in the case 
when A is square and nonsingular. We prove that O(nlogn) samples are sufficient to 
uniquely determine the decomposition with high probability, under the assumption X is 
generated by a Bernoulli-Subgaussian process. 

Our argument for uniqueness suggests a new, efficient dictionary learning algorithm, 
which we call Exact Recovery of Sparsely- Used Dictionaries (ER-SpUD). This algorithm 
solves a sequence of linear programs with varying constraints. We prove that under the 
aforementioned assumptions, the algorithm exactly recovers A and X with high probability. 
This result holds when the expected number of nonzero elements in each column of X is 
at most 0(y/n) and the number of samples p is at least S7(n 2 log 2 ra). To the best of 
our knowledge, this result is the first to demonstrate an efficient algorithm for dictionary 
learning with provable guarantees. 

Moreover, we prove that this result is tight to within a log factor: when the expected 
number of nonzeros in each column is £l(y/n logn), algorithms of this style fail with high 
probability. 

Our algorithm is related to previous proposals by Zibulevsky and Pearlmutter [24] (for 
source separation) and Gottlieb and Neylon [10] (for dictionary learning), but involves 
several new techniques that seem to be important for obtaining provable correct recovery - 
in particular, the use of sample vectors in the constraints. We will describe these differences 
more clearly in Section 5, after introducing our approach. Other related recent proposals 
include [18; 12]. 

The remainder of this paper is organized as follows. In Section 3, we fix our model. 
Section 4 discusses situations in which this problem is well-posed. Building on the intuition 

1. Of course, for some applications, weaker notions than uniqueness may be of interest. For example, 
Vainsencher et. al. [21] give generalization bounds for a learned dictionary A. Compared to the results 
mentioned above, these bounds depend much more gracefully on the dimension and sparsity level. How- 
ever, they do not directly imply that the "true" dictionary A is unique, or that it can be recovered by 
an efficient algorithm. 
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developed in this section, Section 5 introduces the ER-SpUD algorithm for dictionary recov- 
ery. In Section 6, we introduce our main theoretical results, which characterize the regime 
in which ER-SpUD performs correctly. Section 7 describes the key steps in our analysis. 
Technical lemmas and proofs are sketched; for full details please see the full version. Finally, 
in Section 8 we perform experiments corroborating our theory and suggesting the utility of 
our approach. 



2. Notation 

We write for the standard £ p norm of a vector v, and we write ||A<f || for the induced 
operator norm on a matrix M. \\v\\ Q denotes the number of non-zero entries in v. We 
denote the Hadamard (point- wise) product by 0. [n] denotes the first n positive integers, 
{1,2, ... ,n}. For a set of indices /, we let Pj denote the projection matrix onto the subspace 
of vectors supported on indices I, zeroing out the other coordinates. For a matrix X and a 
set of indices J, we let Xj (X J ) denote the submatrix containing just the rows (columns) 
indexed by J. We write the standard basis vector that is non-zero in coordinate i as e^. For 
a matrix X we let row(-X") denote the span of its rows. For a set S, \S\ is its cardinality. 



3. The Probabilistic Models 



We analyze the dictionary learning problem under the assumption that A is an arbitrary 
nonsingular re-by-n matrix, and X is a random sparse n-by-p that follows the following 
probabilistic model: 

Definition 1 We say that X satisfies the Bernoulli-Subgaussian model with parameter 9 
if X = ft R, where ft is an iid Bernoulli(6>) matrix, and R is an independent random 
matrix whose entries are iid symmetric random variables with 



/i = E [ | Rij | ] £ [1/10,1], E[i%] < 1 



and 



' [ \Rij\ >t] < 2exp 



Vt > 0. 



(1) 



(2) 



This model includes a number of special cases of interest - e.g., standard Gaussians and 
Rademachers. The constant 1/10 is not essential to our arguments and is chosen merely 
for convenience. The subgaussian tail inequality (2) implies a number of useful concentra- 
tion properties. In particular, if x\ . . . xn are independent, random variables satisfying an 
inequality of the form (2), then 



N 

1 = 1 



E 



' N 

.i=i 



> t 



<2exp(-^) 



(3) 



We will occasionally refer to the following special case of the Bernoulli-Subgaussian model: 



Definition 2 We say that X satisfies the Bernoulli-Gaussian model with parameter 9 if 
X = ttOR, where ft is an iid Bernoulli(#) matrix, and R is an independent random matrix 
whose entries are iid J\f (0,1). 
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4. When is the Factorization Unique? 

At first glance, it seems the number of samples p required to identify A could be quite 
large. For example, Aharon et. al. view the given data matrix X as having sparse columns, 
each with at most k nonzero entries. If the given samples yj = Axj lie on an arrangement 
of (2) fc-dimensional subspaces range(A/), corresponding to possible support sets /, A is 
identifiable. 

On the other hand, the most immediate lower bound on the number of samples required 
comes from the simple fact that to recover A we need to see at least one linear combination 
involving each of its columns. The "coupon collection" phenomenon tells us that p = 
ft(^ logn) samples are required for this to occur with constant probability, where 9 is the 
probability that an element Xij is nonzero. When 9 is as small as 0(1 /n), this means p 
must be at least proportional to nlogn. Our next result shows that, in fact, this lower 
bound is tight - the problem becomes well-posed once we have observed cnlogn samples. 

Theorem 3 (Uniqueness) Suppose that X = ft R follows the Bernoulli-Subgaussian 
model, and f[Rij = 0] = 0. Then if l/n < 9 < 1/C and p > Cnlogn, with probability at 
least 1 — C'nexp{— c9p} the following holds: 

For any alternative factorization Y = A'X' such that maxj ||e^X'||o < maxj ||e^-X"||o, 
we have A' = All A and X' = A _1 II T X, for some permutation matrix II and 
nonsingular diagonal matrix A. 

Above, c, C, C are absolute constants. 
4.1. Sketch of Proof 

Rather than looking at the problem as one of trying to recover the sparse columns of X, we 
instead try to recover the sparse rows. As X has full row rank with very high probability, 
the following lemma tells us that for any other factorization the row spaces of X, Y and 
X' are likely the same. 

Lemma 4 //rank(X) = n, A is nonsingular, and Y can be decomposed into Y = A'X', 
then the row spaces of X' , X, and Y are the same. 

We will prove that the sparsest vectors in the row-span of Y are the rows of X. As any 
other factorization Y = A'X' will have the same row-span, all of the rows of X' will lie in 
the row-span of Y. This will tell us that they can only be sparse if they are in fact rows of 
X. This is reasonable, since if distinct rows of X have nearly disjoint patterns of nonzeros, 
taking linear combinations of them will increase the number of nonzero entries. 

Lemma 5 Let ft be an n-by-p Bernoulli(#) matrix with l/n < 9 < 1/4. For each set 
SC. [n], let Ts C [p] be the indices of the columns of ft that have at least one non-zero 
entry in some row indexed by S. 

a. For every set S of size 2, 

P[|T 5 | < (4/3)0p] < exp 
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b. For every set S of size a with 3 < a < 1/0 

F[\T S \ < (3o-/8)9p) < exp 

c. For every set S of size a with 1/0 < a, 

F[\T S \ < (l-l/e)p/2] < exp 

Lemma 5 says that every subset of at least two rows of X is likely to be supported on 
many more than Op columns, which is larger than the expected number of nonzeros Op in 
any particular row of X. We show that for any vector a € IR n with support S of size at 
least 2, it is unlikely that a T X is supported on many fewer columns than are in Tg- In the 
next lemma, we call a vector a fully dense if all of its entries are nonzero. 

Lemma 6 For t > 200s ; let fi G {0, 1} SX * be any binary matrix with at least one nonzero 
in each column. Let R G ]R SX * be a random matrix whose entries are iid symmetric random 
variables, with P [Rij = 0] = 0, and let U = ft R. Then, the probability that there exists 
a fully-dense vector a for which \\a T U\\ < t/5 is at most 2~'/ 25 . 

Combining Lemmas 5 and 6, we prove the following. 

Lemma 7 If X = ft R follows the Bernoulli-Subgaussian model, with ¥[Rij = 0] = 0, 
1/n < < 1/C and p > Cnlogn, then the probability that there is a vector a. with support 
of size larger than 1 for which 

||q t x|| < (11/9) 0p 

is at most exp(-cOp). Here, C,c are numerical constants. 

Theorem 3 follows from Lemmas 4 and 7 and the observation that with high probability 
each of the rows of X has at most (10/9) Op nonzeros. We give a formal proof of Theorem 
3 and its supporting lemmas in Appendix A. 

5. Exact Recovery 

Theorem 3 suggests that we can recover X by looking for sparse vectors in the row space of 
Y. Any vector in this space can be generated by taking a linear combination w T Y of the 
rows of Y (here, w T denotes the vector transpose). We arrive at the optimization problem 

minimize ||k; t 1^||o subject to w ^ 0. 

Theorem 3 implies that any solution to this problem must satisfy w T Y = AeJJT for some 
j e [n], A / 0. Unfortunately, both the objective and constraint are nonconvex. We 
therefore replace the 1° norm with its convex envelope, the £ l norm, and prevent w from 
being the zero vector by constraining it to lie in an affine hyperplane {r T w = 1}. This 
gives a linear programming problem of the form 

minimize ||i« t Y]|;l subject to r T w = 1. (4) 

We will prove that this linear program is likely to produce rows of X when we choose r to 
be a column or a sum of two columns of Y. 




1 - l/e)p 
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5.1. Intuition 

To gain more insight into the optimization problem (4), we consider for analysis an equiv- 
alent problem, under the change of variables z = A T w, b = A~ x r: 

minimize ||z T -X"||i subject to b T z = 1. (5) 

When we choose r to be a column of Y, b becomes a column of X. While we do not 
know A or X and so cannot directly solve problem (5), it is equivalent to problem (4): (4) 
recovers a row of X if and only if the solution to (5) is a scaled multiple of a standard basis 
vector: z+ = Xej, for some j, A. 

To get some insight into why this might occur, consider what would happen if X exactly 
preserved the i\ norm: i.e., if = c||z||i for all z for some constant c. The solution 

to (5) would just be the vector z of smallest i 1 norm satisfying b T z = 1, which would 
be ej^/bj^, where j+ is the index of the element of b = A _1 r of largest magnitude. The 
algorithm would simply extract the row of X that is most "preferred" by 6! 

Under the random coefficient models considered here, X approximately preserves the t\ 
norm, but does not exactly preserve it [16]. Our algorithm can tolerate this approximation 
if the largest element of b is significantly larger than the other elements. In this case we 
can still apply the above argument to show that (5) will recover the j*-th row of X. In 
particular, if we let |b|m > 1^1(2) — " " " — l^l(n) De the absolute values of the entries of b 
in decreasing order, we will require both |&|(2)/|&|(i) < 1 — 7 and that the total number of 
nonzeros in b is at most c/9. The gap 7 determines fraction 9 of nonzeros that the algorithm 
can tolerate. 

If the nonzero entries of X are Gaussian, then when we choose r to be a column of Y 
(and thus b = A~ l r to be a column of X), properties of the order statistics of Gaussian 
random vectors imply that our requirements are probably met. In other coefficient models, 
the gap 7 may not be so prominent. For example, if the nonzeros of X are Rademacher (iid 
±1), there is no gap whatsoever between the magnitudes of the largest and second-largest 
elements. For this reason, we instead choose r to be the sum of two columns of Y and thus 
b to be the sum of two columns of X. When 9 < 1/y/n, there is a reasonable chance that 
the support of these two columns overlap in exactly one element, in which case we obtain a 
gap between the magnitudes of the largest two elements in the sum. This modification also 
provides improvements in the Bernoulli-Gaussian model. 

5.2. The Algorithms 

Our algorithms are divided into two stages. In the first stage, we collect many potential 
rows of X by solving problems of the form (4). In the simpler Algorithm ER-SpUD(SC) 
( "single column" ) , we do this by using each column of Y as the constraint vector r in the 
optimization. In the slightly better Algorithm ER-SpUD(DC) ("double column"), we 
pair up all the columns of Y and then substitue the sum of each pair for r. In the second 
stage, we use a greedy algorithm (Algorithm Greedy) to select a subset of n of the rows 
produced. In particular, we choose a linearly independent subset among those with the 
fewest non-zero elements. From the proof of the uniqueness of the decomposition, we know 
with high probability that the rows of X are the sparsest n vectors in row(V). Moreover, 
for p > fi(nlogn), Theorems 8 and 9, along with the coupon collection phenomenon, tell 
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us that a scaled multiple of each of the rows of X is returned by the first phase of our 
algorithm, with high probability. 2 

ER-SpUD(SC): Exact Recovery of Sparsely-Used Dictionaries using single 
columns of Y as constraint vectors. 

For j = 1 . . . p 

Solve miix^ \\w Y\\i subject to (Yej) T w = 1, and set Sj = w T Y. 



ER-SpUD(DC): Exact Recovery of Sparsely-Used Dictionaries using the sum of 
two columns of Y as constraint vectors . 

1. Randomly pair the columns of Y into p/2 groups gj = {Yej 17 Yej 2 }. 

2. For j = 1 . . .p/2 

Let Tj = Yej 1 + Yej 2 , where gj = {Yej 1 , Yej 2 }. 

Solve min^ ||u7 7, y||i subject to rjw = 1, and set Sj = w T Y . 

Greedy: A Greedy Algorithm to Reconstruct X and A. 

1. REQUIRE: S = {s u . . . , s T } C W. 

2. For i = 1 . . . n 

REPEAT 

I argmin S;£l 5 ||s/||oj breaking ties arbitrarily 
Xi = si 
S = S\{s t } 

UNTIL rank([a5i,...,a; i ])=z 

3. Set X = [ Xl ,...,x n ] T , and A = YY T (XY T )- 1 . 

Comparison to Previous Work. The idea of seeking the rows of X sequentially, by 
looking for sparse vectors in row("K), is not new per se. For example, in [24], Zibulevsky 
and Pearlmutter suggested solving a sequence of optimization problems of the form 

minimize Hu^l^lli subject to \\wW2 > 1- 

However, the non-convex constraint in this problem makes it difficult to solve. In more 
recent work, Gottlieb and Neylon [10] suggested using linear constraints as in (4), but 
choosing r from the standard basis vectors e\ . . . e n . 

2. Preconditioning by setting Y p = (YY T )~ 1 ^ 2 Y helps in simulation, while our analysis does not require 
A to be well conditioned. 
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The difference between our algorithm and that of Gottlieb and Neylon — the use of 
columns of the sample matrix Y as linear constraints instead of elementary unit vectors, is 
crucial to the functioning of our algorithm (simulations of their Sparsest Independent Vector 
algorithm are reported below). In fact, there are simple examples of orthonormal matrices A 
for which the algorithm of [10] provably fails, whereas Algorithm ER-SpUD(SC) succeeds 
with high probability. One concrete example of this is a Hadamard matrix: in this case, 
the entries of b = A~ 1 ej all have exactly the same magnitude, and [10] fails because the 
gap between |&|(i) and |6|(2) is zero when r is chosen to be an elementary unit vector. In 
this situation, Algorithm ER-SpUD(DC) still succeeds with high probability. 



6. Main Theoretical Results 

The intuitive explanations in the previous section can be made rigorous. In particular, 
under our random models, we can prove that when the number of samples is reasonably 
large compared to the dimension, (say p ~ n 2 log 2 ra), with high probability in X the 
algorithm will succeed. We conjecture it is possible to decrease the dependency on p to 
0(n log n). 

Theorem 8 (Correct recovery (single-column)) Suppose X is Bernouili(#) — Gaussian. 
Then provided p > c\n 2 log 2 n, and 

2 < e < (6) 

n log n 

with probability at least 1 — Cfp~ 10 , the Algorithm ER-SpUD(SC) recovers all n rows of 
X . That is, all n rows of X are included in the p potential vectors wjY, . . . , w^Y . Above, 
c\, a and Cf are positive numerical constants. 

The upper bound of aj \Jn logn on 9 has two sources: an upper bound of aj^/n is 
imposed by the requirement that b be sparse. An additional factor of ^/log n comes from 
the need for a gap between |b|(i) and |fo|(2) of the k i.i.d. Gaussian random variables. On 
the other hand, using the sum of two columns of Y as r can save the factor of logn in 
the requirement on 9 since the "collision" of non-zero entries in the two columns of X 
creates a larger gap between |&|(i) and |6|(2)- More importantly, the resulting algorithm is 
less dependent on the magnitudes of the nonzero elements in X. The algorithm using a 
single column exploited the fact that there exists a reasonable gap between |6|(i) and |£>|(2), 
whereas the two-column variant ER-SpUD(DC) succeeds even if the nonzeros all have 
the same magnitude. 

Theorem 9 (Correct recovery (two-column)) Suppose X follows the Bernoulli- Subgaussian 
model. Then provided p > c\n 2 log 2 n, and 

2 -<e<^, (7) 

with probability at least 1 — Cfp~ 10 , the Algorithm ER-SpUD(SC) recovers all n rows of 
X . That is, all n rows of X are included in the p potential vectors wjY, . . . , WpY . Above, 
c\, a and Cf are positive numerical constants. 
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Hence, as we choose p to grow faster than n 2 log 2 n, the algorithm will succeed with 
probability approaching one. That the algorithm succeeds is interesting, perhaps even 
unexpected. There is potentially a great deal of symmetry in the problem - all of the 
rows of X might have similar ^ 1 -norm. The vectors r break this symmetry, preferring one 
particular solution at each step, at least in the regime where X is sparse. To be precise, 
the expected number of nonzero entries in each column must be bounded by y/n logn. 

It is natural to wonder whether this is an artifact of the analysis, or whether such a bound 
is necessary. We can prove that for Algorithm ER-SpUD(SC), the sparsity demands in 
Theorem 9 cannot be improved by more than a factor of y/\og n. Consider the optimization 
problem (5). One can show that for each i, ||e^JT||i ~ 9p. Hence, if we set z = e^/bj^, 
where j* is the index of the largest element of b in magnitude, then 

||„T Y |. II^^H 1 _ r B P 

||o||oo \/logn 

If we consider the alternative solution v = sign(b)/||6||i, a calculation shows that 

||u T X||i « C'p/y/n. 

Hence, if 9 > cydogn/n for sufficiently large c, the second solution will have smaller 
objective function. These calculations are carried through rigorously in the full version, 
giving: 

Theorem 10 If X follows the Bernoulli- Subgaussian model with 



jPlogn 



n 



with (3 > /3q, and the number of samples p > cnlogn, then the probability that solving the 
optimization problem 

minimize ||io T y||i subject to r T w = 1 (8) 
with r = Yej recovers one of the rows of X is at most 

exp (— cp) + exp ^— 3/3 v /n l°g n^j + 4exp [—dOp + logn) + 2n 1 ™ 5/3 (9) 

above, Po,c,c' are positive numerical constants. 

This implies that the result in Theorem 8 is nearly the best possible for this algorithm, 
at least in terms of its demands on 9. A nearly identical result can be proved with r = 
Yei + Yej the sum of two rows, implying that similar limitations apply to the two-column 
version of the algorithm. 
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7. Sketch of the Analysis 

In this section, we sketch the arguments used to prove Theorem 8. The proof of Theorem 
9 is similar. The arguments for both of these results are carried through rigorously in 
Appendix B. At a high level, our argument follows the intuition of Section 5, using the 
order statistics and the sparsity property of b to argue that the solution must recover a row 
of X . We say that a vector is A;-sparse if it has at most k non-zero entries. Our goal is to 
show that z± is 1-sparse. We find it convenient to do this in two steps. 

We first argue that the solution z+ to (5) must be supported on indices that are non-zero 
in b, so z is at least as sparse as b, say -^/n-sparse in our case. Using this result, we restrict 
our attention to a submatrix of y/n rows of X, and prove that for this restricted problem, 
when the gap 1 — |b|(2)/|fr|(i) is large enough, the solution z± is in fact 1-sparse, and we 
recover a row of X. 

Proof solution is sparse. We first show that with high probability, the solution z± to 
(5) is supported only on the non-zero indices in b. Let J denote the indices of the s non-zero 
entries of \b\, and let S = {j \ Xjj ^ 0} C [p], i.e., the indices of the nonzero columns 
in Xj. Let zq = PjZ* be the restriction of z+ to those coordinates indexed by J, and 
z\ = z+ — zq. By definition, zq is supported on J and z\ on J c . Moreover, zq is feasible 
for Problem (5). We will show that it has at least as low an objective function value as z+, 
and thus conclude that Z\ must be zero. Write 

||z^jc||i = ||,z^jr s '||i + ||i > ||z^x s '||i — n^f x^ii! + ll^i 1 .^ ||i 

= \z\X\x - 2\\z\X s \\ x + \z\X\x, (10) 

where we have used the triangle inequality and the fact that ZqX s = 0. In expectation 
we have that 

\\z^X\\t > \\zlX\\ 1 + {p-2\S\)E[\\zJX\\ 1 ] > \\z%X\\ 1 + c(p-2\S\)x/efii\\z 1 \\ 1 , (11) 

where the last inequality requires On > 2. 

So as long as p — 2\S\ > 0, zq has lower expected objective value. To prove that this 
happens with high probability, we first upper bound IS*! by the number of nonzeros in Xj, 
which in expectation is Osp. As long as p — 2(1 + 5)6sp = p(l — c'Os) > 0, or equivalently 
s < c s /6 for some constant c s , we have [|^JC[|i > ||z^X||i. In the following lemma, we 
make this argument formal by proving concentration around the expectation. 

Lemma 11 Suppose that X satisfies the Bernoulli-Subgaussian model. There exists a nu- 
merical constant C > 0, such that if On > 2 and 

p > Cn 2 log 2 n (12) 

then with probability at least 1 — 3p~ 10 , the random matrix X has the following property: 

(PI) For every b satisfying 
problem 

minimize 
has supp(z*) C supp(b). 



b\\o < 1/80, any solution z* to the optimization 
z T X\\ 1 subject to b T z = 1 (13) 
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Note in problem (5), b = A~ 1 r. If we choose r = Yei, then b = A~ 1 Ye{ = Xei, and 
E[||6||o] = On. A Chernoff bound then tells us that with high probability z* is supported 
on no more than 29n entries, i.e., s < 28n. Thus as long as 29n < C2/O, i.e., 9 < ce/y/n, we 
have || 2* ||o < 2#n = cgy/n. 



The solution in Xj: If we restrict our attention to the induced s-by-p submatrix Xj, 
we observe that X j is incredibly sparse - most of the columns have at most one nonzero 
entry. Arguing as we did in the first step, let j* denote the index of the largest entry of |6j|, 
and let S = {j \ X 7^ 0} C [p], i.e., the indices of the nonzero entries in the j*-th 
row of Xj. Without loss of generality, let's assume bj* = 1. For any z, write zq = Pj*z 
and z\ = z — zq. Clearly zq is supported on the j*-th entry and Z\ on the rest. As in the 
first step, 

||2 T Xj||i > - 2\\z'(X S j\\ 1 + H^XjHl (14) 

By restricting our attention to 1-sparse columns of Xj, we prove that with high probability 

\\zJXj\\ l >fi9p(l-se)(l-e) 2 \\z 1 \\ 1 . 
We prove that with high probability the second term of (14) satisfies 

\\zJXj, s \\i < (l + e)/i0 2 p||zi||i. 

For the first term, we show 

Wz^Xj^ > \\ejtXj\l! - l&^ZilUXjIli > HeJ^XjIl! - |^zi|(l + e)fi9p. 

If |6|( 2 )/|6|(i) < 1 - 7, then \b T jZ i\ < (1 - 7)lki||l- 

In Lemma 12, we combine these inequalities to show that if if 9 < cy/s, then 

||(2 + z 1 ) T X J \\ 1 > ||eJ*A:j||i + ^p(l-7)||zi||i. (15) 

Since ej* is a feasible solution to Problem 5 with a lower objective value as long as z\ ^ 0, 
we know ej* is the only optimal solution. The following lemma makes this precise. 

Lemma 12 Suppose that X follows the Bernoulli-Subgaussian model. There exist positive 
numerical constants cx,C2 such that the following holds. For any 7 > and s G Z + such 
that 9s < 7/8 and p is sufficiently large: 

c is log n ) p > C2_ 

9~f 2 ' J logp — #7 2 ' 



then with probability at least 1 — 4p 10 , the random matrix X has the following property: 

(P2) For every J E (^) and every b£l s satisfying |&|(2)/|&|(i) — 1 — 7> ^ e 
solution to the restricted problem, 

minimize \\z T X j^\\i subject to b T z = 1, (16) 

is unique, 1-sparse, and is supported on the index of the largest entry ofb. 

Once we know that a column of Y provides us with a constant probability of recovering 
one row of X, we know that we need only use 0(n log n) columns to recover all the rows 
of X with high probability. We give proofs of the above lemmas in Section B. Section B.3 
shows how to put them together to prove Theorem 8. 
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8. Simulations 



ER-SpUD(proj): Exact Recovery of Sparsely-Used Dictionaries with Iterative 
Projections . 

S <- {0} C R n . 

For i = 1 . . . n 

For j = 1 . . . p 

Find W{j G argmin^ [|wj*Y'[|i subject to (Yej) T P s ±w = 1. 
Wi <(— argmin l(;=lUai ... )lu . T ||w*"K||o, breaking ties arbitrarily. 
5^—5© span(iOj). 

X <- W. 

a «- yr(ir)^ 1 

In this section we systematically evaluate our algorithm, and compare it with the state- 
of-the-art dictionary learning algorithms, including K-SVD [1], online dictionary learning 
[15], SIV [10], and the relative Newton method for source separation [23]. The first two 
methods are not limited to square dictionaries, while the final two methods, like ours, exploit 
properties of the square case. The method of [23] is similar in provenance to the incremental 
nonconvex approach of [24], but seeks to recover all of the rows of X simultaneously, by 
seeking a local minimum of a larger nonconvex problem. We found in the experiments that 
a slight variant of the greedy ER-SPUD algorithm, we call the ER-SPUD(proj), works even 
better than the greedy scheme. 3 And thus we also add its result to the comparison list. As 
our emphasis in this paper is mostly on correctness of the solution, we modify the default 
settings of these packages to obtain more accurate results (and hence a fairer comparison). 
For K-SVD, we use high accuracy mode, and switch the number of iterations from 10 to 
30. Similarly, for relative Newton, we allow 1,000 iterations. For online dictionary learning, 
we allow 1,000. We observed diminishing returns beyond these numbers. Since K-SVD and 
online dictionary learning tend to get stuck at local optimum, for each trial we restart K- 
SVD and Online learning algorithm 5 times with randomized initializations and report the 
best performance. We measure accuracy in terms of the relative error, after permutation- 
scale ambiguity has been removed: 

fe(A,A) = min ||AAn — 
n,A 

Phase transition graph. In our experiments we have chosen A to be a an n-by-n matrix 
of independent Gaussian random variables. The coefficient matrix X is n-by-p, where p = 
5n log e n. Each column of X has k randomly chosen non-zero entries. In our experiments 
we have varied n between 10 and 60 and k between 1 and 10. Figure 1 shows the results for 
each method, with the average relative error reported in greyscale. White means zero error 

3. Again, preconditioning by setting Y p = (YY T )~ 1 /' 2 Y helps in simulation. 
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(a) ER-SpUD(SC) (b) ER-SpUD(proj) (c) SIV 




(d) K-SVD (e) Online (/) Rel. Newton 



Figure 1: Mean relative errors over 10 trials, with varying support k (y-axis, increase from bottom 
to top) and basis size n(x-axis, increase from left to right). Here, p = 5nlog e n. Our 
algorithm using a column of Y as r (ER-SpUD(SC)), our algorithm with iterative pro- 
jections (ER-SpUD(proj)), SIV [10], K-SVD [1], online dictionary learning [15], and the 
relative Newton method for source separation [23]. 



and black is 1. The best performing algorithm is ER-SpUD with iterative projections, which 
solves almost all the cases except when n = 10 and k > 6. For the other algorithm, When 
n is small, the relative Newton method appears to be able to handle a denser X, while as n 
grows large, the greedy ER-SpUD is more precise. In fact, empirically the phase transition 
between success and failure for ER-SpUD is quite sharp - problems below the boundary are 
solved to high numerical accuracy, while beyond the boundary the algorithm breaks down. 
In contrast, both online dictionary learning and relative Newton exhibit neither the same 
accuracy, nor the same sharp transition to failure - even in the black region of the graph, 
they still return solutions that are not completely wrong. The breakdown boundary of K- 
SVD is clear compared to online learning and relative Newton. As an active set algorithm, 
when it reaches a correct solution, the numerical accuracy is quite high. However, in our 
simulations we observe that both K-SVD and online learning may be trapped into a local 
optimum even for relatively sparse problems. 

9. Discussion 

The main contribution of this work is a dictionary learning algorithm with provable per- 
formance guarantees under a random coefficient model. To our knowledge, this result is 
the first of its kind. However, it has two clear limitations: the algorithm requires that 
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the reconstruction be exact, i.e., Y = AX and it requires A to be square. It would be 
interesting to address both of these issues (see also [2] for investigation in this direction). 
Finally, while our results pertain to a specific coefficient model, our analysis generalizes to 
other distributions. Seeking meaningful, deterministic assumptions on X that will allow 
correct recovery is another interesting direction for future work. 
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Appendix A. Proof of Uniqueness 

In this section we prove our upper bound on the number of samples for which the decom- 
position of Y into AX with sparse X is unique up to scaling and permutation. We begin 
by recording the proofs of several lemmas from Section 4. 
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A.l. Proof of Lemma 4 

Proof Since rank(X) = n, we know rank(A') > rank(l^) = rank(A) = n. Since both A 
and A' are nonsingular, the row spaces of X' and X are the same as that of Y. ■ 



A. 2. Proof of Lemma 5 

Proof First consider sets S of two rows. The expected number of columns that have 
non-zero entries in at least one of these two rows is 

p(l - (1 - Of) = p(26 - 6 2 ) > (3/2)p0, 

for 6 < 1/2. Part a now follows from a Chernoff bound. 

For part b, is a > 3 and aO < 1, we observe that for every S 



e\t s \ = p-{i-ey P > (oe-&)e 2 )p = (i-^e^aep > ^ 



where the inequalities follow from aQ < 1. Part b now follows from a Chernoff bound. 
For part c, if a9 > 1, for every S of size a we have 

E \T S \ > (1 - e- a9 )p > (1 - e^p. 

As before, the result follows from a Chernoff bound. 



Definition 13 (fully dense vector) We call a vector a € K" fully dense if for all i G 
[n], cti / 0. 

A. 3. Proof of Lemma 7 

We use the following theorem of Erdos. 

Theorem 14 ([6]) For every k > 2 and nonzero real numbers z%, ... ,Zk, 



Ziri = 



< 2" 



k 

[k/2\ 



< 1/2, 



where each ri is chosen independently from ±1, 



Lemma 15 For b > s, let H G M sxb 6e any matrix with at least one nonzero in each 
column. Let R be an s-by-b matrix with Rademacher random entries, and let U = H 5] . 
Then, the probability that the left nullspace of U contains a fully dense vector is at most 



-b+s \og(e 2 b/s) 
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Proof As in the preceding lemma, we let U = [u±\ . . . denote the columns of U and 
for each j 6 [b], we let Nj be the left nullspace of [«i| . . . \v,j]. We will show that it is very 
unlikely that Ni, contains a fully dense vector. 

To this end, we show that if iVj_i contains a fully dense vector, then with probability at 
least 1/2 the dimension of Nj is less than the dimension of Nj—±. To be concrete, assume 
that the first j — 1 columns of X! have been fixed and that iVj_x contains a fully dense 
vector. Let a be any such vector. If Uj contains only one non-zero entry, then cx T Uj ^ 
and so the dimension of Nj is less than the dimension of iV,-_i. If Uj contains more than 
one non-zero entry, each of its non-zero entries are random Rademacher random variables. 
So, Theorem 14 implies that the probability over the choice of entries in the j'th column of 
SI that a T Uj = is at most one-half. So, with probability at least 1/2 the dimension of Nj 
is less than the dimension of Nj—i. 

To finish the proof, we observe that the dimension of the nullspaces cannot decrease 
more than s times. In particular, for iV& to contain a fully dense vector, there must be at 
least b — s columns for which the dimension of the nullspace does not decrease. Let F C [b] 
have size b — s. The probability that for each j S F that Nj—i contains a fully dense vector 
and that the dimension of Nj equals the dimension of Nj—i is at most 2~ b+s ~ 1 . Taking a 
union bound over the choices for F, we see that the probability that iVj contains a fully 
dense vector is at most 

^ ^2~ 6+s = f^2~ 6+s < f—~\ 2~ b+s < 2~ 6+s + sl °g( efe / s ) < 2~ b+slo g( e2b / s ^ 



Proof [Proof of Lemma 6] Notice that if X = Q, R is Bernoulli-Subgaussian, then 
because the entries of R are symmetric random variables, X is equal in distribution to 
ft R X, where S is an independent iid Rademacher matrix. We will apply Lemma 15 
with H = ft Q R. 

If there is a fully-dense vector a for which ||o! T £/|| < i/5, then there is a subset of at 
least b = At/5 columns of U for which a is in the nullspace of the restriction of U to those 
columns. By Lemma 15, the probability that this happens for any particular subset of b 
columns is at most 

2~ b+s\oge 2 b/s ^ 2~ 4i/5+ s log(e 2 t/s) 

Taking a union bound over the subsets of b columns, we see that the probability that this 
can happen is at most 

( * | 2 ~ 4i/5+sloge2 * / ' s < 2 a722 '2~'( 4 / 5- ( s /*) log(e2 '/ s ^ < 2*( a722 ~ a8+0 - 0365 ) < 9~*/ 25 

where in the first inequality we bound the binomial coefficient using the exponential of the 
corresponding binary entropy function, and in the second inequality we exploit s/t < 1/200. 
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A. 4. Proof of Lemma 7 

Proof Rather than considering vectors, we will consider the sets on which they are sup- 
ported. So, let S C [n] and let a = \S\. We first consider the case when 17 < a < 1/0. 
Let T be the set of columns of X that have non-zero entries in the rows indexed by S. Let 
t = \T\. By Lemma 5, 

P [t < (3/8)<r0p] < exp(-a0p/64). 

Given that t > (3/8)a9p, Lemma 6 tells us that the probability that there is a vector a 
with support exactly S for which 

||q t X|| < (ll/9)0p < (3/4O)a0p 

is at most 

exp(-(3/2OO)cr0p). 

Taking a union bound over all sets S of size a, we see that the probability that there vector 
a of support size a such that ||q; t X|| < (ll/9)0p is at most 

Q (exp(-(3/2OO)cr0p) + exp(-cr0p/64)) < exp(-ca0p), 

for some constant c given that p > Cn log n for a sufficiently large C. 

For a > 1/0, we may follow a similar argument to show that the probability that there 
is a vector cc with support size a for which ||a r X|| < (ll/9)0p is at most 

exp(-cp), 

for some other constant c. Summing these bounds over all a between 17 and n, we see that 
the probability that there exists a vector a with support of size at least 17 such that such 
that ||a T X|| < (ll/9)0p is at most 

exp(— c9p), 

for some constant c. 

To finish, we sketch a proof of how we handle the sets of support between 2 and 17. For 
a this small and for sufficiently small relative to a (that is smaller than some constant 
depending on a), each of the columns in T probably has exactly one non-zero entry. Again 
applying a Chernoff bound and a union bound over the choices of S, we can show that 
with probability 1 — exp(— c9p) for every vector a with support of size between 2 and 17, 
||a T X|| > (5/4)0p. ■ 



A. 5. Proof of Theorem 3 

Proof From Lemma 7 we know that with probability at most exp(— c9p), any dense linear 
combination of two or more rows of X has at least (ll/9)0p nonzeros. Hence, the n rows 
of X are the sparsest directions in the row space of Y. 
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A Chernoff bound shows that the probability that any row of X has more than 

(lO/9)0p 

non-zero entries is at most 




Hence, with the stated probability, the rows of X are the n sparsest vectors in row(X). 

On the aforementioned event of probability at least 1 — exp(— c9p), X has no left null 
vectors with more than one nonzero entry. So, as long as all of the rows of X are nonzero, 
X will have no nonzero vectors in its left nullspace. With probability at least 1 — n(l — 6) p > 
1 — nexp(— cp), all of the rows of X are nonzero, and so row(X) = iow(Y) = row(X'). 

This, together with our previous observations implies that every vector in row(X') is a 
scalar multiple of a row of X , from which uniqueness follows. Summing failure probabilities 
gives the quoted bound. ■ 



Appendix B. Proof of Correct Recovery 

Our analysis will proceed under the following probabilistic assumption on the coefficient 
matrix X: 

Notation. Below, we will let ||Af|| r i = max, ||e^7Vf||i, where the are the standard 
basis vectors. That is to say, || • || r i is the maximum row I 1 norm. This is equal to the 
I 1 — > i 1 operator norm of M T . In particular, for all v, M, \\v T M\\i < \\v\\i\\M\\ r i. 

B.l. Proof of Lemma 11 

Proof We will invoke a technical lemma (Lemma 17) which applies to Bernoulli-Subgaussian 
matrices whose elements are bounded almost surely. For this, we define a truncation oper- 
ator T T : M. nx P R nx P via 

Mij \Mij\ < t 
else 

That is, T T simply sets to zero all elements that are larger than r in magnitude. We will 
choose r = ^/24\ogp and set 

X' = T r [X] = ftQT T [R] = ftQR'. (17) 

The elements of R' are iid symmetric random variables. They are bounded by r almost 
surely, and have variance at most 1. Moreover, 

/■oo poo rco 

/i' = ^[\R'ij\] = / F[\R! i:j \>t)dt = / F[\Rij\ > t] dt - P[\Rij\>t]dt 

Jt = Jt = Jt=T 

> ^~ 2 J t expT-Mdt > n-2p- 12 > 1/20. 



(T T [M])ij = 
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The final bound follows from provided the constant C in the statement of the lemma is 
sufficiently large. Of course, since \x < 1, we also have fjf < 1. 

The random matrix X' is equal to X with very high probability. Let 

£ x = event{X = X'}. (18) 

We have that 

F[£ c x ] = P[3(i,j) | |Xy|>r] < 2npexp(-y) . (19) 

Ensuring that C is a large constant (say, C > 2 suffices), we have p > n. Since r = v / 241og~p, 
P[££] < 2p- 10 . 

For each J C [n] of size at most s = 1/89, we introduce two "good" events, £s(I) and 
En (I)- We will show that on £x n£s(I) n£/v(I), for any b supported on /, and any optimal 
solution z*, supp(z*) C I. Hence, the desired property will hold for all sparse b on 

£ g ood = £x n fl £<?(/) n £ N {i). (20) 
|/|<« 

For fixed I, write T{I) = {j \ X itj = V % G I}, and set S(I) = [p] \ T(I). That is 
to say, T is the set of indices of columns of X whose support is contained in I c , and S is 
its complement (indices of those columns that have a nonzero somewhere in I) . The event 
£s(I) will be the event that S(I) is not too large: 

£ S (I) = event{\S(I)\ <p/4}. (21) 

The event £n(I) will be one on which the following holds: 

V«eR n - |!| , ||« T X^||i-2||t> T A:^ 5(J) ||i > cip/i'^Hi. 
Since // > 1/20, this implies 

Vv£M. n -W, \\v T X' IC J 1 -2\\v T X' ICyS{I) \\ 1 > cap ^ Hi. (22) 

Obviously, on fx n£jv(7), the same bound holds with X' replaced by X . Lemma 17 shows 
that provided S is not too large, £n{I) is likely to occur: F[£n(I) \ £s{I)] is large. 

On £x n £s(I) fl £n(I), we have inequality (22). We show that this implies that if b is 
sparse, for any solution z+, supp(z^) C supp(6). Consider any b with supp(fo) C I, and any 
putative solution z+ to the optimization problem (5). If supp(2 >r ) C /, we are done. If not, 
let z £ R n such that 

M - { |2 o'- * • ' 23 > 

and set zi = — zq. Notice that since b 1 Z\ = 0, zq is also feasible for (5). We prove that 
under the stated hypotheses z\ = 0. 



Form a matrix X' s G M nxp via 
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and set X' T = X' — X' s . We use the following two facts: First, since S and T are disjoint, 
for any vector q, \\q\\i = ||<Zs||i + ||QtIIi- Second, by construction of T, z^Xt = 0. Hence, 
we can bound the objective function at z* = zo + z\ below, as 

IIOso + *i) T *lli = 11(^0 + z 1 ) T X s \\i + 11(^0 + zxfXrWt 

> \\z^X s \\i - \\zJXs\\i + II*i-*t||i. (25) 

Since ZqX = Zq(X$ + Xt) = ZqX$, this bound is equivalent to 

|| (z + *i) T X||i > Wz^XWx + \\z\ X\\i - 1\\z[ X s ||i (26) 

Noting that z\ is supported on 7 C , (22) implies that if z\ 7^ 0, ||zo"X||i is strictly smaller 
than HzJ'Xlli = IK^o + zi) T X||i. Hence, zo is a feasible solution with objective strictly 
smaller than that of z*, contradicting optimality of z*. To complete the proof, we will show 
that P[£good] is large. 

Probability. The subset S is a random variable, which depends only on the rows of Xj^ 
of X indexed by /. For any fixed /, |<S| ~ Binomial(p, A), with A = 1 — (1 — We have 
< A < \I\6, where the upper bound uses convexity of (1 — 9)^. From our assumption on 
s, A < \I\6 < 1/8. Applying a Chernoff bound, we have 

F[\S\>p/4] < ¥[\S\>2Xp] < exp(-^f) < exp [-% J . (27) 



3 J ~ 1 V 3 , 

Hence, with probability at least 1 — exp ( — t)> we nave l^l < 

Since X' is iid and S depends only on X'j ^, conditioned on S, X\ c ^ is still iid Bernoulli- 
Subgaussian. Applying Lemma 17 to X' JC conditioned on S gives 

F[£ N (I) C \£ S (I)] < expf + « log (Cn ^ogp)) ■ (28) 

In this bound, we have used that fjf G [1/20, 1], r = y / 241ogp, and the fact that the bound 
in Lemma 17 is monotonically increasing in n to simplify the failure probability. Moreover, 
we have 

F[£ s (I)n£ N (I)] = 1-F[£ N (I) C \£ S (I)]F[£ S (I)] -F[£ S (I) C ] 

> l-F[£ N (iy\£ s (I)} - F[£s(I) c } 

> 1 - exp (-— ^j== + nl °g {C n V^Zpj^j ~ exp f~y V 
Let s max = [§5 J denote the largest value of |J| allowed by the conditions of the lemma. 

F[q ood ] < F[£ c x ]+ £ F[(£ s (i)n£ N (l)) c ] 

Ic[n], \I\<s max 

< m\ + — ( {»? (-|) + - P (--^= + n log (Cn^H)) } 
2 J -- 1 " + S „„(jJ {«*(-§;) +exp(--|= + „log(c„ N ^,))} 



< 



< 2p~ 10 + exp( % — + C"n log (n log p) ) 

V nVlogp / 
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for appropriate constant C" . Under the conditions of the lemma, provided C is large enough, 
the final term can be bounded by p~ 10 , giving the result. ■ 



Lemma 16 Suppose x = fi R G W 1 with fi iid Bernoulli(#) ; R an independent random 
vector with iid symmetric entries, n9 > 2, and [i = E[ \Rij\ ] < +oo. Then for all v G R n ; 



E f|v T aj|l > —\l -\\v\\ . 
u 4 V n 



(29) 



Proof Let z = v/\\v\U, and write 



■[IE' 



E 



[\^2 n i R j z i\ 



71 

> inf EflV a-ii.C,- ll = inf EflV^&C, 

||<|| 1= i Vj^ J C>o,EiO=i L "V 

where the final equality holds due to the fact that sign(£j)iL is equal to -Rj in distribution. 
Notice that E | YJ ■ f2ji?j£j| is a convex function of ((j), and is invariant to permutations. 



Hence, this function is minimized at the point d = C2 



Cn = 1/n, and 



£ > n _1 E[|J^%i2j| 



[IE 1 



(30) 



Let T = #{j I £lj = 1}. For fixed O, Y^- fij-Rj is a sum of symmetric random variables. 
Thus 



(31) 



where (e) is an independent sequence of Rademacher (iid ±1) random variables. By the 
Khintchine inequality, 



E^l^^li?^ 



> -^E R \\n R\\ 2 > -^=E R \\fl QR\\i = 1L\[T]2. (32) 
V 2 v2T 



Hence 



ri Z — ^ 



(33) 



s=0 



Notice that T ~ Binomial(n, 9), and hence, with probability at least 1/2, T > \n0\ > nO/2, 
where the last inequality holds due to the assumption n9 > 2. Plugging in to (33), we 
obtain that 



e > 



4 V n 



(34) 
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Lemma 17 Suppose that X follows the Bernoulli- Subgaussian model, and further that for 
each < r almost surely. Let X$ be a column submatrix indexed by a fixed 

set S C [p] of size \S\ < |. There exist positive numerical constants c, C such that the 
following: 



\v T X\\ 1 -2\\v T X s \\ 1 > 



HP 



32 V n 

holds simultaneously for all v £ W n , on an event 8^ with 



(35) 



< exp 



+ n log ( C max { — , 1 



n(l + r/i) 

Proof Consider a fixed vector z of unit i 1 norm. Note that 



nr 



|z^.X^||i — 2||z'^-X" i 5'||i 



(36) 



ies c its 

is a sum of p independent random variables. Each has absolute value bounded by 

„T 



and second moment 



z Xi\ < z i LX" oo < t 



\w{\z T Xi\) < E[(z T Xi) 2 } < 9\\z\\\ < 9. 



There are p such random variables, so the overall variance is bounded by p9. The expectation 
of is 

E [\\z T X\\x - 2||z T X 5 || 1 ] = (p - 2\S\)E\z T X 1 \ =(p- 2\S\)fi > !f, 

where fl = E|z T Xi|. We apply Bernstein's inequality to bound the deviation below the 
expectation: 



|z T X||i-2||z T X 5 ||i < y~* 



< exp 



29p + 2rt/3 



(37) 



We will set t 



16 V n 



p, and notice that by Lemma 16, fi > 



4 V n 



Hence, we obtain 



\z T X\\ 1 -2\\z T X s \\ 1 < j^\ -P 

16 V n 



< exp 



cifi 2 9p 2 /n 



29p + c 2 pnTy / 9/n 



(38) 



Using that 9/n < # 2 /2 (from the assumption 9n > 2), [i < 1, and simplifying, we can obtain 
a more appealing form: 



\z T X\\ 1 -2\\z T X s \\i < t^\I-V 

16 V n 



< exp 



cfj, 2 p 



n(l + r/i) / ' 



(39) 
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This bound holds for any fixed vector of unit i l norm. We need a bound that holds 
simultaneously for all such z. For this, we employ a discretization argument. This argument 
requires two probabilistic ingredients: bounds for each z in a particular net, and a bound 
on the overall I 1 operator norm of X T , which allows us to move from a discrete set of z to 
the set of all z. We provide the second ingredient first: 

Now, let N be an e-net for the unit "1-sphere" r = {x \ \\x\\i = 1}. That is, for all 
Ijccj^ = 1, there is a z 6 N such that ||x — z|L < e. For example, we could take 

N = {y/ \n/e] :yeZ n ,\\y\\ 1 = \n/e]}. 

With foresight, we choose e = ^-W- > 4tt^-. The inequality follows from 9n>2. Using 
standard arguments (see, e.g. [20, Page 15]), one can show 

m < 2 «^ n / £ 1 + n " X \ < (2e(\l/e] + l)) n . 



So, 



lo.o |;V | < //log ( C'max^ ^,1 J> j . (10) 



for appropriate constant C. 
For each z € N, set 



£ z = event I ||z T X||i - 2||z T X s ||i > ^ ( ' ' ) 



: A\ < <■*!>(- J - (12) 



Our previous efforts show that for each z, 

/?( 1 + r//) 
Set 

£n= (~)£ z . (43) 

On £jv, consider any uGT, choose z £ N with — z||i < e, and set A = — z, then 
\\v T X\\ 1 -2\\v T X s \\ 1 = \\(z + AfXsc\\i-\\(z + A) T X s \\i 

> ||^ T ^5 c ||l — ||2 T ^s||l — ||A T X5c||i — ||A T X5||i 

= ||z T X||i - 2||z T X 5 ||i - ||A T X'||i 

> ||z T X||i - 2||z T X 5 ||i - ||A||i||X T ||i 

u [e 

16 ~ £T \!f) ]/n P 

&\fip. (44) 
32 V n 
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This bound holds all v G T (above, we have used that for v £T, \\v\\i = 1). By homogeneity, 
whenever this inequality holds over T, it holds over all of M. n we have 



||u r X||i-2||u T X s ||i > £v-pHi. ( 45 ) 

62 \ n 

To complete the proof, note that the failure probability is bounded as 

< E p [^] ^ ^ |xexp (-n(iT^)) 



( c/j 2 p ( fnr 

< exp ; + n log L max < — , 1 

V n(l + r/u) V I 



B.2. Proof of Lemma 12 

Proof For each j G [n], set 1} = | = 0}. Set 

Ojj = | X u £ 0, and = 0, Vj' G J \ {j}} . 
Consider the following conditions: 

\\X\\ rl <(l + e)iJ0p, (46) 
v i G N, ll- X '[n]\{j},T/l|rl < a^P, (47) 

VJef^VieJ, ||X/,n.Jll > /3/x0p. (48) 

We show that when these three conditions hold for appropriate e,a, f3 > 0, the property 
(P2) will be satisfied, and any solution to a restricted subproblem will be 1-sparse. Indeed, 
fix J of size s and nonzero b G M. n , with supp(6) C J. Let j* G J denote the index of the 
largest element of b in magnitude. 

Consider any v whose support is contained in J, and which satisfies b T v = 1. Then we 
can write v = vq + v\, with supp(^o) Q {j*} and supp(i>i) C J \ {j*}. We have 

||« T Xj||i = Wv^Xj + vJXjW, 

> WvZXjWx + \\v\ X JjT ||i - ||^X J)T c t Id. (49) 

J 3 

Above, we have used that vq is supported only on j*, supp(i^X) = IJ*, and applied the 
triangle inequality. 

First without loss of generality, we assume by normalization bj* = 1. We know that 
vq = (1 — b T vi)ej*, and the vector ej* is well-defined and feasible. We will show that in 
fact, this vector is optimal. Indeed, from (49), we have 

H^XjHx > \\(l-b T v 1 )ef,Xj\\ 1 + \\v^Xj^\\ 1 -\\vfX J>T cJ 1 

> \\eJ*Xj\\i — \\X jWriWbj^j^WooWvxWi + WvJXj^Tj* Id ~~ ||' y l||l||-X'j,Tf lt llrl 

> ||eJ*Xj|d + \\vJX JtT ^ Id " (ll^jlln x (1 - 7 ) + || r i) ||«i|d 
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In the last simplification, we have used that |&|(2)/H(i) < 1 — 7- We can lower bound 
Wv^X^Tj* ||i as follows: for each i G J, consider those columns indexed by j G Oj,i. For 
such j, 

suppOi) n supp(Xjj) = {i}, 

and so 

ll«i'-Xj,nj, 4 ||i = biWIII^n^lli- 

Notice that i ^ i' , tljj n fij^/ = 0, and that for j G J, and i / j, tljj C Tj. So, finally, we 
obtain 

K^>||i > ^I^WHlX^Jl!. (50) 

ieJ 

Plugging in the above bounds, we have that 

||u T X||i > \\ej ir X\\ 1 + (fi + 'y-a-E-l)\\v 1 \\ifi0p. (51) 
Hence, provided /3 — a — e > 1 — 7, Bj* achieves a strictly smaller objective function than v. 

Probability. We will show that the desired events hold with high probability, with (3 = 
1 — 7/4, and a = e = 7/4. Using Lemma 18, (46) holds with probability at least 

( cj 2 9p , 
1 — 4 exp I — h log n 

For the second condition (47), fix any j. Notice that for any (j' , i) the events Xj> i 7^ and 
i G Tj are independent, and have probability 6 2 . Moreover, for fixed j', the collection of 
all such events (for varying column i) is mutually independent. Therefore, for each j E [n], 
|| -X"[ji]\{j},Tf [|rl is equal in distribution to the rl norm of an (n — 1) x p matrix whose rows 
are iid # 2 -Bernoulli-Subgaussian. 4 Hence, by Lemma 18, we have 

l X H\{i},r;l|r-i > (l + 5)n9 2 p\ < 4 exp (-c5 2 8 2 p + log(n - 1)) . (52) 

To realize our choice of a, we need 9(1 + 5) = 7/4; we therefore set 5 = 9~ l (7/4 — 9) > 
9~ 1 ^y/8. Plugging in, and taking a union bound over j shows that (47) holds with probability 
at least 

l-4exp(-^^ + 21ogn) . (53) 

Finally, consider (48). Fix j and J. Notice for i G [p], the events i G are independent, 
and occur with probability 9' = 9(1 — 6) s ~ 1 > 9(1 — 6s). Hence H-X^n^ . ||i is distributed as 
the i 1 norm of a 1 x p iid ^'-Bernoulli-Gaussian vector. Again using Lemma 18, we have 

IP [\\X jfljij ||i < (1 - 5)i29'p] < 4 exp (-c5 2 9'p) . (54) 

To achieve our desired bound, we require (1 — 5)9' > (39; a sufficient condition for this is 
(1 — 5)(1 — 9s) > 1 — 7/4, or equivalently, 

5(1 -9s) < 7/4 -9s. 

4. Distinct rows are not independent, since they depend on common events i G Tj, but this will cause no 
problem in the argument. 
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Using the bound 9s < 7/8, we find that it suffices to set 5 = 7/8. We therefore obtain 

P [\\Xj,njJl < Ptfp] < 4exp • (55) 

Finally, using 9' > 9(1 - 9s) > 9/2, we obtain 

P [ll^njll < MP] < 4exp ("^) • (56) 

Taking a union bound over the s( n ) pairs j,J and summing failure probabilities, we see 
that the desired property holds on the complement of an event of probability at most 

. (n\ ( cj 2 9p\ A ( cy 2 9p \ / c^ 2 9p 
4s I J exp I — I + 4 exp I — + 2 log n j + 4 exp I — + log n 

< / A»\ 2 A ( H 2 ' 

< 4 s + rT + n x exp ' 



sV 7 1 V 128 



Using that log (™) < slogn and 9 > 2/n, it is not difficult to show that under our hypothe- 
ses, the failure probability is bounded by Ap~ 10 . ■ 



Lemma 18 Let X be an n x p random matrix, such that the marginal distributions of the 
rows Xi follow the Bernoulli- Subgaussian model. Then for any < 5 < 1 we have 

(1 -<&)// 0p< \\X\\ r i < (l + 5)fi9p (57) 

with probability at least 

1 — 4 exp iy—c5 2 9p + log n) . 
where c is a positive numerical constant. 

Proof Let Xi = Ri denote the i-th row of X. We have ||-X"||ri = maxj g j n ] ||Xj||i. 
For any fixed i, \\Xi\\i = ^ \Xij\ = Yjj \^ijR%j\, and E[||Xj||i] = Y^j^\ x ij\ = V®P> where 
[i = E,[\Rij\] is as specified in the definition of the Bernoulli-Subgaussian model. 

Let Ti = = 1}. Since Tj ~ Binomial(p, 9), for < 5i < 1, the Chernoff bound 

yields 

P [Ti > (1 + ^Ap] < exp (-^) , (58) 

P [Ti < (1 - .JOflp] < exp (-^) • (59) 

Since the Rij are subgaussian, so is \Rij\. Conditioned on fij, ||-X"i||i = X^jeTj * s a sum 
of independent subgaussian random variables. By (2), 

F[\\\X i \\ 1 -iiT i \>5 2 »e P \n i ] < 2 exp f- y j . (60) 
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Whenever € Op x [1 — 5%, 1 + 6i], the conditional probability is bounded above by 



exp 



" 2(1-5!) 



So, unconditionally, 

£ /x0p x [1 - Si - S 2 , 1 + <5i + 5 2 ) 
c\b\iiOp 



< 2 exp 



[I -Si) 



+ exp 



5ft 



+ exp 



5lt 



Set 5i = = 5/2, the fact that /x is bounded below by a constant, combine exponential 
terms, and take a union bound over the n rows to complete the proof. ■ 



B.3. Proof of Correct Recovery Theorem 8 

Proof If n = 1, the result is immediate. Suppose n > 2. We will invoke Lemmas 11 and 12. 
The conditions of Lemma 11 are satisfied immediately. In Lemma 12, choose 7 = j^^, with 
/3 smaller than the numerical constant ao m the statement of Lemma 24. Take s = \68ri\ . 
Since On > 2, s < 70n. The condition 0s < 7/8 is satisfied as long as 

0' < JL = P , 

56n 56nlogn 

This is satisfied provided the numerical constant a in the statement of the theorem is 
sufficiently small. Finally, it is easy to check that p satisfies the requirements of Lemma 12, 
as long as the constant c\ in the statement of the theorem is sufficiently large. 

So, with probability at least 1 — 7p -10 , the matrix X satisfies properties (P1)-(P2) 
defined in Lemmas 11 and 12, respectively. Consider the optimization problem 

minimize Hu^l^Hi subject to r T w = 1, (61) 

with r = Yej the j-th column of Y. This problem recovers the i-th row of X if the solution 
is unique, and A T w ir is supported only on entry i. This occurs if and only if the solution 
z* to the modified problem 

minimize ||z T X||i subject to b T z = 1 (62) 

with b = A~ 1 r = Xej is unique and supported only on entry i. 

Provided the matrix X satsifies properties (P1)-(P2), solving problem (61) with r = 
Yej recovers some row of X whenever (i) 1 < ||Jfej||o < s and (ii) |fr|(2)/H(i) < 1 — 7- Let 

Si{j) = event{ HXejIlo > }, (63) 
£ 2 (j) = event{ ||X ej || < s }, (64) 
£ 3 (j) = event{ |6| (2) /|6| (1) < 1 - 7 }• (65) 

Let £(i,j) be the event that these three properties are satisfied, and the largest entry of 
b = Xej occurs in the i-th entry: 

£{i,j) = £x(j) n £ 2 (j) n £ 3 (j) n event { |6| (1) = \b t \ } . 
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If the matrix X satisfies (P1)-(P2), then on £(i,j), (61) recovers the i-th row of X. 
Moreover, because £(i,j) only depends on the j-th column of X, the events £ (i, 1) . . . £(i,p) 
are mutually independent. By symmetry, for each i, 

F[£(i,j)] = -F[s 1 {j)ns 2 (j)ne 3 (y)]. 

n 

The random variable ||-X"ej||o is distributed as a Binomial(n, 9). So, P[£i(j) c ] = (1 — 9) n < 
(1 — 2/n) n < e~ 2 . The binomial random variable ||.X"ej[|o has expectation 9n. Since s > 69n, 
by the Markov inequality, F[£ 2 (j) c ] < 1/6. Finally, by Lemma 24, F[£ 3 (j) c \ £{\ < 1/2. 
Moreover, 

P[£iC?)n5 2 tf)n£ 3 0')] > i - P^i(i) c ] - P[^(i) c ] - F[£ 3 (j) c \ £{\ 

> 1 - e~ 2 - 1/6 - 1/2 = C- 
The constant £ is larger than zero (one can estimate £ ~ .198). For each 

F[£(i,j)] > -• 
n 

Hence, the probability that we fail to recover all n rows of X is bounded by 

n 

F[X does not satisfy (P1)-(P2)] + 5^P[Djf (i, jf\ < 7p" 10 + n(l - C/^) p 



i=l 



< 7p 10 + exp ( — — + log n 

n 



Provided ^ > logn + 10 log p, the exponential term is bounded by p 10 . When c\ is chosen 
to be a sufficiently large numerical constant, this is satisfied. ■ 



B.4. The Two-Column Case 

The proof of Theorem 9 follows along very similar lines to Theorem 8. The main difference 
is in the analysis of the gaps. 

Proof We will apply Lemmas 11 and 12. For Lemma 12, we will set 7 = 1/2, and 
s = Y19n + 1. Then 9s = Y19 2 n. Under our assumptions, provided the numerical constant 
a is small enough, 9s < 7/8. Moreover, the hypotheses of Lemma 11 are satisfied, and so 
X has property (PI) with probability at least 1 — 3p~ 10 . Provided a is sufficiently small, 
s < 1/89, as demanded by Lemma 11. 

For each j S \p\, let Qj = supp(Xej). Let 

£n,j,k(i) = event U Ofe| < s, and Qj flflt = {i}} . 

Hence, on £n,j,k(i), the vectors Xej and Xe^ overlap only in entry i. The next event will be 
the event that the i-th entries of Xej and Xe^ are the two largest entries in the combined 
vector 

" Xe< 



h 



3 k 



Xe k 
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and their signs agree. More formally: 

£x,j,k{i) = event {{\h jk \ {1) , \h jk \ (2) } = {X i:j ,X ik } and sign(Xjj) = sign(X ifc ) } . 

If we set r = Yej + Ye k , and hence b = A~ 1 r = Xej + Xe k , then on £x,j,k{i) H £nj, k (i), 
the largest entry of b occurs at index i, and |&|(2)/|&|(i) < 1/2. 

Hence, if X satisfies (P1)-(P2), on £n,j,k(i) ^£x,j,k(i), the optimization 

minimize ||tu T "F||i subject to r T w = 1 (66) 

with r = Yej + Ye k recovers the i-th row of X. The overall probability that we fail to 
recover some row of X is bounded by 

n 

F[X does not satisfy (P1)-(P2)] + J> [nf/ 2 1 (fnj(i) 1 fc(i)( i ) n ^xj(i) 1 fc(i)( i )) t 

i=i 

Let O'- = fij \ and fi^ = fife \ {z}. Then we have 



[ f nj,fc(») 



> 



[jenjandieOfe] P [ Ufi' fc |<s-1 and ^ n fif k 

s-1 



1 



W > ^ 



|n»l> 



[n;nn' fc = 0] -: 
> e 2 ((i - 2 ) n_1 - 1/3) , 

where in the final line we have used the Markov inequality: 

2E[|n; 



-p[n;ni2' fc ^0] 
ItVJ > ^ 



in'-l > — 



< 



s-1 



20(n-l) 1 
s-1 ~ 6' 



Since 9 < a/y/n, we have (1 — 2 ) n 1 > 1 — a 2 . Provided a is sufficiently small, this quantity 
is at least 2/3. Hence, we obtain 

P[&Ufc(«)] > ^ 

For £x,j,k{i)-> we calculate 
and 



> 



2(£) - 2(12fln + l) 2 ' 



F[%,t(!)n^(i)] > 



for some numerical constant £ > 0. Hence, the overall probability of failure is at most 



7p- 10 + n I 1 



p/2 



< 7p~ 10 + exp 



2n 2 



+ logn 



(67) 

(68) 
(69) 

(70) 



Provided J^- > logn + 10 log this quantity can be bounded by p 10 . Under our hypothe- 
ses, this is satisfied. ■ 



30 



ER-SPUD 



Appendix C. Upper Bounds: Proof of Theorem 10 

Proof For technical reasons, it will be convenient to assume that the entries of the random 
matrix X are bounded by r almost surely. We will set X' = T T X, and set r = logp, so 
that with probability at least 1 — 2p~ 10 , X' = X. We then prove the result for X' . Notice 
that we may write X' = fi R! , with R' iid subgaussian, and fjf = E[|i^-|] > 1/20. 

Applying a change of variables z = A T w, b = A~ 1 r, we may analyze the equivalent 
optimization problem 



minimize \\z X'\\\ subject to b z = 1. 



(71) 



Let v = nmfp- This vector is feasible for (71). We will compare the objective function 
obtained at z = v to that obtained at any feasible one-sparse vector z = e-i/bi. Note that 



|sign(6) T X'|| 1 = 



k=l 



X 'ik s[ S n (bi 



i=l 



— Y^Qk — Qj + Qk- 

k k^j 



The random variable Qj is just the &■ norm of X'p which is bounded by nr. The random 
variables (Qk)k^j are conditionally independent given b. We have 



E [Ql 



011611 



By Bernstein's inequality, 







+ t\X> 


_k+j 







< exp 



-t 2 



(72) 



^2^||X;.||o + 2r||X;.|| /3 / 
If the constant c in the statement of the theorem is sufficiently large, then pO > c\t. We 



also have E [Q^ \ X'j\ < */6\\Xj\\o. Simplifying and setting t = p*/ 0\\Xj\\o, we obtain 



J^Qk >2p^||x;.|| 1 x) 

k^j 



< exp (— c 2 p) ■ 



(73) 



; n 



The variable Qj is bounded by nr. Moreover, a Chernoff bound shows that 

P[||Xj-|| > 40n] < exp(-3<9n) < exp (-3Py/n\og; 
So, with overall probability at least 1 — exp (— C2p) + exp (— 3/3^n logn), 

||X /T sign(5)||i < AOpy/n + nr < 59py/n, 



(74) 



(75) 



where in the final inequality, we have used our lower bound on p. Using that b = X 'ej is 
Bernoulli-Subgaussian, we next apply Lemma 18 with 5 = 1/2 to show that 



||6||i < \y!dn\ < 4exp(— c^On) < 4exp(— c^^/ flnlogn) 



(76) 
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On the other hand, by Lemma 18, 

P min||efX'||i < 

i 

and using the subgaussian tail bound, 



< 4exp(— C30p + logn) 



P[||6||oo >*] < 2nexp(-t 2 /2). 

Set t = -y/lOlogn to get that ||6||oo < \/10/31ogn with probability at least 2n 1-5 ^. 
Hence, with overall probability at least 



1 — exp (— C2p) — exp ( — 3/3 yn logn ) — 4exp (—c^Op + logn) — 2n 



,1-5/3 



we have 



while 



/T sign(6) 



X 



,T<2i 



X 



II \-' T „ II 

\\X e< i 



§9p^/n csp 



fi9n 



n 



9 HP 



> 



c^9p 



> 



n 



■n 



(77) 



(78) 



(79) 



(80) 



2V10Vlog^ 

If (3 is sufficiently large, (80) is larger than (79), and the algorithm does not recover any 
row of X' . Since with probability at least 1 — 2p~ 10 , X 1 = X, the same holds for X, with 
an additional failure probability of 2p~ 10 . ■ 



Appendix D. Gaps in Gaussian Random Vectors 



In this section, we consider a (i-dimensional random vector r, with entries iid Af(0, 1). We 
let 

> s(2) > > s(d) 

denote the order statistics of \r\. We will make heavy use of the following facts about 
Gaussian random variables, which may be found in [7, Section VII. 1]. 

Lemma 19 Let x be a Gaussian random variable with mean and variance 1. Then, for 
every t > 



p(t) > F [x > t] > 



1 1 

t~1? 



where 



p(t) 



'2tt 



exp (-t 2 /2) . 



Lemma 20 For any d > 2, 



s(l)>4 V / l^d 



< d 



-3 
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Proof Set t = 4\/log d in Lemma 19 to show that 

P[|r(i)|>*]< J— d-\ (81) 
4v2vrlogd 

The denominator is larger than one for any d > 2; bounding it by 1, and taking a union 
bound over the d elements of r gives the result. ■ 



Lemma 21 For d larger than some constant 



sd) < ^ 



< exp 



4 v / log d 



Proof Set 

t = -^logd 

For d sufficiently large, we may use Lemma 19 to show that the probability that a Gaussian 
random variable of variance 1 and any mean has absolute value greater than t is at least 



4 v / log d 

Thus, the probability that every entry of r has absolute value less than t is at most 

d 7 / 8 \ 



< exp - 



4^1og d 



We now examine the gap between the largest and second-largest entry of r. We require 
the following fundamental fact about Gaussian random variables. 

Lemma 22 Let x be a Gaussian random variable with variance 1 and arbitrary mean. 
Then for every t > and a > 0, 

P [\x\ <t + a\\x\>t] < 3amax(t,3). 

Proof Assume without loss of generality that the mean of x is positive. Then, 

P [t < \x\ < t + a] 



[\x\ < t + a\ \x\ > t] 



< 



F[\x\ > t] 
2P [t < x < t + a] 

Ppj > t] 
2P [t < x < t + a] 



~ P[x>t] 
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One can show that this ratio is maximized when the mean of x is in fact 0, given that it is 
non-negative. Similarly, the ratio is monotone increasing with t. So, if t < 3, we will upper 
bound this probability by the bound we obtain when t = 3. When the mean of x is and 
t > 3, Lemma 19 tells us that 

P [x > t] > —L^- ( 1 - i ] exp(-i 2 /2) > — exp(-t 2 /2). 



We then also have 

F[t<x<t + a]= / exp(-x 2 /2)dx < exp(-i 2 /2). 

The lemma now follows from 2-9/8 < 3. 



Lemma 23 Let d>2. Then for every a > 0, 

P [a(l) - s(2) < a] < 12a\/W + d~ 2 . 

Proof Let Mj be the event that is the largest entry of q in absolute value. As the 
events Mj are disjoint, we know that the sum of their probabilities is 1. Let Ri be the event 
that max.^ \q{j)\ < 4\/\nd. From Lemma 20 we know that the probability of not(i?j) is 
at most 1/d 3 . Let Gj be the event that Mj holds but that the gap between q(i) and the 
second-largest entry in absolute value is at most a. From Lemma 22, we know that 

P [Gi\Ri and Mj] < UaVhui. 

The lemma now follows from the following computation. 

P[*(l)- S (2)<a]=J>[Gy 

i 

= P [Gi and Mj] 

i 

= ^P [Gi and Mj and Ri) + ^P [G { and Mj and not(i?j)] . 

i i 

We have 



and 



^P [d and Mj and not(i?,)] < J^P [ not(i?j)] < 1/d 2 , 

i i 

^2 P [°i alld M i and = P t Mi aIld P and ^ 

i i 

< F i M i\ p [Gt\Mi and Jfc] 

i 

< ^P[Mj] 12aVhui 



12avW. 
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Lemma 24 There exists ao > such that for any a < ao and 2 < d < n the following 
holds. If r is a d- dimensional random vector with independent standard Gaussian entries 

s(l) > s(2) > > s(d) 

are the order statistics of \r\, then 



1 _s{2) < a 



si 



log n 



Proof For any a > 0, t > 0, we have 



i _s(2) > a 



logn 



> P 

> 1 

> 1 



«(1) -«(2) > 



< 



at 



1 



(82) 



and s(l) < t\f\og 



n 



8(1) -8(2) < 
8(1) -8(2) < 



at 



at 



-y/Iogn 

Setting t = 4 and applying Lemmas 23 and 20, we have 



s(l) > tv/bg" 

8(1) > ty/lojd 



1 _s_(2) > a 



s(l) logn 



> 1 - 48au 



'logfi 
logn 



eT 2 - 



Choosing ao sufficiently small (and noting that for d > 2, d 2 + d 3 < 3/8) completes the 
proof. ■ 
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